---
title: "BigFig2"
output: html_document
---


#Larvae Standard DESEQ analysis using karyotype 
```{r}
library(DESeq2)
library("magrittr")
library("dplyr")
library("ggplot2")
library("reshape2")
library("RColorBrewer")
library("gplots")
library("ggbio")

matrix = read.table("larvae.counts.matrix")
samples = read.table("larvae.txt", sep="\t", header=TRUE, row.names = 1)
samples$ID = rownames(samples)

# remove all rows with less than 10 total counts
larvae_new = matrix %>% mutate(Total = dplyr::select(., L1:YBB2) %>% rowSums(na.rm = TRUE))
rownames(larvae_new) = rownames(matrix)
larvae_new = subset(larvae_new, larvae_new[ , 29] > 10, drop=FALSE ) 
larvae_new = larvae_new[,1:28]

matrix = as.matrix(larvae_new)
matrix2 = round(matrix)


dds <- DESeqDataSetFromMatrix(countData = matrix2,
                              colData = samples,
                              design = ~ Genotype + Population)
#remove L3v1

dds <- dds[,-11]

#remove SKBB1
dds <- dds[,-24]
pure = dds[,c(1,6,7,11,18:26)]
cross = dds[,c(2:5,8:10,12:17)]



dds = DESeq(dds)

resultsNames(dds)
plotDispEsts(dds)

```





#Larvae look at overall sample relationships using matrix and pca
```{r}
vsd <- vst(dds, blind=FALSE)
rld <- rlog(dds, blind=FALSE)


new = assay(vsd)

greeks = c(expression(paste(alpha,alpha)),expression(paste(beta,beta)), expression(paste(alpha,beta)))

ggplotColours <- function(n = 6, h = c(0, 360) + 15){
  if ((diff(h) %% 360) < 1) h[2] <- h[2] - 360/n
  hcl(h = (seq(h[1], h[2], length = n)), c = 100, l = 65)
}

library(ggfortify)


vsd$Karyotype = vsd$Genotype

pcaData <- plotPCA(vsd, intgroup=c("Karyotype"), returnData=TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))
pcaData$Karyotype = ordered(pcaData$Karyotype, levels=c("AA", "BB", "AB"))

goodL = ggplot(pcaData, aes(PC1, PC2, color=Karyotype)) +
  geom_point(size=6) +
  xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  ylab(paste0("PC2: ",percentVar[2],"% variance")) + scale_colour_manual(values=ggplotColours(n=3),labels=greeks) +
  coord_fixed()  + theme_bw() + theme(legend.position = "none", axis.title.x = element_text(size=14), axis.title.y = element_text(size=14), axis.text.x = element_text(size=12), axis.text.y = element_text(size=12)) 


pdf("pca_larvae.pdf",width=6,height=6)
goodL
dev.off()

##Create a fake larvae sex for making the legend

vsd2 = vsd
vsd2$Sex = c(rep("Male",13),rep("Female",13))


pcaData <- plotPCA(vsd2, intgroup=c("Karyotype", "Sex"), returnData=TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))
pcaData$Karyotype = ordered(pcaData$Karyotype, levels=c("AA", "BB", "AB"))

legend = ggplot(pcaData, aes(PC1, PC2, color=Karyotype, shape=Sex)) +
  geom_point(size=4) +
  xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  ylab(paste0("PC2: ",percentVar[2],"% variance")) + scale_colour_manual(values=ggplotColours(n=3),labels=greeks) +
  coord_fixed()


```

#adults
```{r setup, include=FALSE}
library(DESeq2)
library("magrittr")
library("dplyr")
library("ggplot2")
library("reshape2")
library("RColorBrewer")
library("gplots")

matrix = read.table("adult2.counts.matrix")
samples = read.table("adults2.txt", sep="\t", header=TRUE)
samples$ID = samples$good


adults_new = matrix %>% mutate(Total = dplyr::select(., P10F:P9FD) %>% rowSums(na.rm = TRUE))
rownames(adults_new) = rownames(matrix)
adults_new = subset(adults_new, adults_new[ , 18] > 10, drop=FALSE ) 
adults_new = adults_new[,1:17]


matrix = as.matrix(adults_new)
matrix2 = round(matrix)

```

#run DESEQ for all adults
```{r}
dds <- DESeqDataSetFromMatrix(countData = matrix2,
                              colData = samples,
                              design = ~ Sex + Genotype + Sex*Genotype)
dds = DESeq(dds)
resultsNames(dds)
plotDispEsts(dds)
```

#Adults look at overall sample relationships using matrix and pca
```{r}
vsd <- vst(dds, blind=FALSE)


greeks = c(expression(paste(alpha,alpha)), expression(paste(beta,beta)))

ggplotColours <- function(n = 6, h = c(0, 360) + 15){
  if ((diff(h) %% 360) < 1) h[2] <- h[2] - 360/n
  hcl(h = (seq(h[1], h[2], length = n)), c = 100, l = 65)
}

vsd$Karyotype = vsd$Genotype

pcaData <- plotPCA(vsd, intgroup=c("Karyotype", "Sex"), returnData=TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))
goodA = ggplot(pcaData, aes(PC1, PC2, color=Karyotype, shape=Sex)) +
  geom_point(size=6) +
  xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  ylab(paste0("PC2: ",percentVar[2],"% variance")) + scale_colour_manual(values=ggplotColours(n=3),labels=greeks) +
  coord_fixed() + theme_bw()+ theme(legend.position = "none",axis.title.x = element_text(size=14), axis.title.y = element_text(size=14), axis.text.x = element_text(size=12), axis.text.y = element_text(size=12)) 




pdf("pca_adults.pdf",width=6,height=6)
goodA
dev.off()



```


```{r}
library(cowplot)
library(ggpubr)
library(patchwork)

fig1 = png::readPNG('Figure 1.png') 

im_A <- ggplot() + background_image(fig1) + 
    theme(plot.margin = margin(t=0, l=0, r=0, b=0, unit = "cm"))

legend_b2 <- get_legend(legend +theme(legend.direction="horizontal",legend.position="bottom"))



top = ggarrange(im_A, goodL, ncol=2, labels=c("A","B"),widths = c(1, 0.9))
bottom = plot_grid(goodA, legend_b2, ncol =1, rel_heights  = c(1, .15),align = "hv", labels=("C") )

P2 = ggarrange(top,bottom, ncol=1, nrow=2)


pdf("Figure_1.pdf",width=12,height=12)
P2
dev.off()


```

